﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using TomanuExtensions;
using System.Diagnostics;
using RaytracerLib.MathLib.RootFinders;

namespace RaytracerLib.MathLib.Polynomials
{
    public struct PolynomialSolverAprox
    {
        private readonly SturmSequence m_sturm_sequence;
        private RootFinder m_root_finder;

        public PolynomialSolverAprox(Polynomial a_polynomial, RootFinder a_root_finder)
        {
            m_sturm_sequence = new SturmSequence(a_polynomial);
            m_root_finder = a_root_finder;
        }

        public double FindOneRootInRange(double a_a, double a_b, int a_num_roots_at_a)
        {
            for (int i = 0; i < Constants.STURM_FIND_ONE_ROOT_AREA_MAXIMUM_ITERATIONS; i++)
            {
                double mid = (a_a + a_b) / 2;
                int atmid = m_sturm_sequence.NumberOfSignChangesAt(mid);
                int nroots = a_num_roots_at_a - atmid;

                if (nroots == 1)
                    return m_root_finder.FindRoot(a_a, mid, m_sturm_sequence[0]);
                else if (nroots == 0)
                    a_a = mid;
                else
                    a_b = mid;
            }

            return a_a;
        }

        public double SolveOne()
        {
            int nroots_at_a = m_sturm_sequence.NumberOfSignChangesAt(Constants.MINIMAL_DISTANT);
            int nroots_at_pos_inf = m_sturm_sequence.NumberOfSignChangesAtPositiveInfinity();

            int nroots = nroots_at_a - nroots_at_pos_inf;

            if (nroots == 0)
                return Double.PositiveInfinity;

            double a = Constants.MINIMAL_DISTANT;
            double b = a + 1.5;
            int nroots_at_b = nroots_at_a;

            for (int its = 0; its < Constants.STURM_ANY_ROOTS_MAXIMUM_ITERATIONS; its++)
            {
                nroots_at_b = m_sturm_sequence.NumberOfSignChangesAt(b);

                if (nroots_at_a - nroots_at_b >= 1)
                    break;

                b *= 2.0;
            }

            if (nroots_at_a == nroots_at_b)
                return Double.PositiveInfinity;

            return FindOneRootInRange(a, b, nroots_at_a);
        }
    }
}
